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Abstract. The technological innovation of the experimental techniques involved in 
molecular biology allows nowadays to perform experiments almost unimaginable a few 
decades ago. Among these a particularly interesting role is played by single molecule 
manipulation experiments, aiming at understanding the mechanical properties of the 
fundamental bricks of life. DNA stretching experiments show for instance that both 
single stranded and double stranded DNA (respectively ss-DNA and ds-DNA) have first 
an elastic elongation, due to entropic effects, follow by an enthalpic regime in which the 
molecule's constituents themselves start to elongate. Moreover, when the stretching 
force reaches values around 65 pN ds-DNA can suffer a second type of transition, after 
which the molecule reaches an elongation of about 1.7 times its contour length. 
Besides those fascinating experiments it is necessary to formulate mathematical models 
to give a quantitative flavour to the already large amount of experimental results. 
These models need to be realistic enough to describe the complexity observed in 
the experiments yet sufficiently simple to be conducive to fitting the experimental 
data. So far all the models proposed to reproduce experimental measurements of DNA 
elongation failed to satisfy both grounds, either lacking some fundamental properties 
of the molecule or being too involved to analyse experimental data. 
In this paper we reconsider a model proposed a few years ago by Storm and Nelson 
[15j which represented an interesting step towards a theory that could synthesise the 
properties of accuracy and simplicity. We solve the model using the cavity method, 
avoiding any mathematical approximation, yielding results that not only fit successfully 
experimental measurements performed on single and double stranded DNA but also 
produce sensible physical values of the parameters of the model. Furthermore, to 
confirm the mathematical quality of the approach, our findings are also compared to 
Monte Carlo simulations. 
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1. Introduction 

The impressive development of experimental techniques in molecular biology made 
possible a brand new class of experiments, among which a particularly important role is 
played by the ones of single molecule manipulation. Such experiments aim to understand 
the response of the very fundamental bricks of living beings to a mechanical stimulus 
(e.g. DNA pulling experiments can give some insights on the mechanical behaviour of 
the Nucleic Acid during replication). 

The pioneering works of Bustamante [1] in the early 90s opened a research field that 
uncovered many peculiar features of biological molecules under tension, showing, for 
example, that DNA and RNA have an elastic response to small pulling forces and a 
resistance to bending [21 [3] , implying that the correlation along the molecular chain is 
non negligible. 

Some of these experiments showed as well another interesting property of double 
stranded DNA (dsDNA): it has been seen that when B-DNA (i.e. the DNA in its 
well-known helical configuration) is stretched by forces of around 65 pN, it suffers a 
structural change and suddenly the molecule reaches an extension of about 1.7 times 
its natural length [U |5]. The nature of this overstretched DNA (which has been named 
S-DNA) stimulated a long debate in the last years [6l [7] , as two different visions of the 
hydrogen bonds breaking could explain the transition [HI [HI HO]. Indeed, while it was 
clear that the double helix unwinds and the base pair align with the stretching force, 
it was not fully understood whether the B-S transition was caused by a complete melt 
into bubbles of the molecule (strand unpeeling) or rather a local bond breaking (a pure 
transition to a new state), until a coherent explanation was found [11]. In this last work 
it is highlighted how both processes can happen according to experimental conditions, 
which can lead to either a fast transition to S-DNA or a single stranded molecule. 
In order to understand better the mechanical properties of bio-molecules, it is essential 
to complement the experimental approach with theoretical models which are able to 
reproduce the main features of those systems. The simplest model that captures 
at least some basic mechanical features is the so-called Freely Jointed Chain (FJC) 
[I2] . This model depicts the molecule as a chain of non interacting units which is 
stretched by an external force and where the monomers are discrete segments with 
an orientation in space. FJC model is qualitatively in agreement with the typical 
experimental measurements (e.g. in the regime of small forces), with some differences 
arising from the fact that pulled DNA extends even more than its contour length. In 
the intermediate regime of forces the FJC solution is however quantitatively different 
from the real elongation curves, showing then that the intermolecular interaction plays 
an important role. 

A second model that incorporate more features of real polymers is the Worm Like Chain 
(WLC) [13], the standard choice to fit the experimental measurements. In this model, 
the molecule is seen as a whole continuous body, where the interaction is felt as a 
bending stiffness, which characterises the correlation at zero force. The WLC solution 
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is a good improvement respect to the FJC, but an accurate comparison of the WLC 
result with real data still shows discrepancies between the model and actual systems 
[2]. One main difference comes from the fact that the fundamental building blocks of 
DNA behave extensibly, which results for high force into an elongated DNA longer than 
its native length at rest. This feature is then normally taken into account considering 
the building blocks as springs and introducing a so called "stretching modulus" [H]. 
However this care is not enough, since the WLC can be a good approximation as long 
as the correlation length is much longer than the base pair length in DNA. 
Some years ago Storm and Nelson [15] introduced the Discrete Persistent Chain (DPC) 
model, where both the coarse-grained nature (like in the FJC) and the bending stiffness 
(as in the WLC) of real molecules are incorporated. In this model the DNA pulling 
experiments are described using a Heisenberg-like Hamiltonian, where the stretching 
force is an external field and the interaction along the chain makes the polymer stiff. 
Their solution (where again a stretching modulus was considered) was in fair agreement 
with the experimental data of [I6] and lead to a fit which quantitatively improved the 
previous WLC results. Unfortunately, the author's nice modelling is handicapped by an 
analytical treatment relying in a series of approximations which yielded unreasonable 
values as, for instance, the monomer's length being much smaller than the real base pair 
size for ss-DNA. 

It is important to note that all these are phenomenological models and to improve upon 
them one needs firstly to understand what they are actually modelling: what does the 
discrete segment capture from the real DNA molecule? Does it model, for instance, a 
base pair or a bigger part of the DNA? In order to understand this model we revisit 
here the DPC model and fit it with experimental measurements. 

2. The Model 

2.1. Single stranded DNA 

In the DPC model [151 [13 [H] the DNA is considered to be a chain of interacting 
monomers with Hamiltonian 

N N-1 

- /3H{i) = fbB ^U-z + JbY^ U ■ U+i (1) 

1 = 1 4 = 1 

where t = (ti, . . .t^) , f is the stretching force, bs the monomers' length and Jb the 
interaction between the monomers, which tunes the bending stiffness of the chain. As we 
expect the segments to be aligned and to resist bending, Jb is chosen as a ferromagnetic 
coupling. Vectors z and ti lie on the unit sphere, the former pointing to the stretching 
direction, while the latter are aligned with the direction in space of the monomers. Note 
that the Hamiltonian ([1]) is simply a Heisenberg model [19] in external field fbB- The 
free energy per monomer F is 

-/3F = llog l^^d^te-^^W, (2) 
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where the subscript denotes the integral over the unit sphere. As the fundamental 
blocs of real DNA are extensible themselves, the stretched DNA is longer than its contour 
length. To capture this feature we introduce a "stretching modulus" and rewrite the 
monomer's length as Bb = b^^^ {1 + f /Y) with Y the Young modulus. The chain's 
elongation then reads 

i(/;^)^^^|i:('-.-> (3) 

i=l 

with cb = b^B (1 + 2//y), where (■ ■ ■) is the usual thermal average, and fi = [b^B ^ Jb, Y) 
are the parameters of the model. Note that the dependence of the elongation on the 
force is explicitly through cb and implicitly through the thermal average. 



2.2. Double stranded DNA 

The preceding model can be extended so as to capture the behaviour of double stranded 
B-DNA when pulled by strong forces. Since the first experiments by Smith et al |1] 
and Cluzel et al [5] it has been observed that double stranded B-DNA undergoes a 
very sharp transition when stretched by a force around 65 pN. When the pulling force 
reaches this critical value, the molecule extends its contour length of about 1.7 times 
in a very abrupt manner, showing a highly cooperational transition. The state of this 
over-stretched DNA is called S-DNA and its emergence has been tightly related to 
the experimental setup [H]. To capture this we reconsider the Hamiltonian ([T]) and, 
as in [15], we allow the monomers to be into two possible states: native B-DNA and 
denaturated S-DNA [181 CSl EO] • This yields the following Hamiltonian: 

N N 

- {t, cr) = / ^ b^Ii ■ z + -fB^S^,,B + Yl i-^^r^^^+i^i ■ ^Vi + ea„<x,+i) (4) 

i=l i=l j=l 

where cXi G {i?, S} represents the state of monomer i. Here b^ represents the monomer's 
length in state a and J„y[= J^'^a) the interaction between two neighbouring monomers 
in states a and a'. In what follows we assume e^a = and denote J„ = We also 
assume, as before, that the monomer's length of the native state is 6^ = b^^ (1 + f IY\ 
while the one in the denaturated state is rigid. The chain's elongation reads: 

N h ^ 

i=l i=l 

where the parameters of this second model are ^ = {fo^'', Y, bs, Jb, Jbs, Js, 1b, ^bs}- 



3. Parameter fitting 

The parameters /x can be fit to a set of experimental measurements {{f!l''^^\ La^^^)}'^=i 
by using Newton's method, viz. 

^(^^^) = f,i^)-H.V^x\ i = 0,l,... (6) 
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Here H and V^x^ represent the Hessian and gradient of the cost function 



X 



a=l 



respectively. They read as follows: 



5' 



(7) 



(8) 



a=l 



a=l 



z,j = 1,2,3. (9) 



As it is expected that lI''''^'' ^ L (^fa"^^^] Atj, the second term in the Hessian is usually 
discarded and the latter is approximated by 



(10) 



a=l 



Note that the expressions for the gradient V^x^ ^iiid the approximated Hessian 
depend on the derivatives of the elongation with respect to the parameters of the model, 
and these are simply correlation functions. For instance, in model ([T]), derivatives of L 
with respect to Jb and 6^ (derivatives with respect to and Y follow from hs via the 
chain rule), yield the following correlation functions: 



dL 



N N-\ 



B 



N 



dL _ fcB_ 



N 



N 



5^((tV5) (4 ■5))^, (11) 



1=1 k=l i,k=l 

where (■ ■ ■)^ denotes connected correlation in the thermal average. 
Thus in order to minimise efficiently the cost function by Newton's method, we 
need to find a way to deal with the thermal average appearing in eqs. f l5|TT]) and the 
corresponding expressions for model (jlj). Luckily the two models introduced are simple 
enough that, as we will see, a simple set of equations can be found to perform thermal 
averages efficiently. 



4. Analytical solution with cavity method 

4.1. Single stranded DNA 

Due to the one-dimensional nature of the problem, the solution of the DPC model is 
usually sought via the transfer matrix method, possibly expanding the eigenfunctions in 
spherical harmonics. However, while in the Heisenberg model with zero field the transfer 
matrix can be easily diagonalised with such an expansion [21] , in this case the two terms 
in ([1]) do not commute and hence a variational method is generally adopted. Another 
source of approximations is to consider the asymptotic behaviour of the elongation 
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L{f; /x) for small and high forces and use the subsequent expressions (see Appendix A ) 
to fit the data. As we want to fit the model in the whole range of forces to assess 
its validity, avoiding, in turn, unnecessary approximations, we proceed as follows: we 
use the cavity method [22] to write exact solutions for the constrained local partition 
function and the chain's elongation; then we evaluate the correlation functions (ITTil 
involved in the minimisation of the cost function ([7]) by writing down equations for the 
propagation of a perturbation in the cavity equations. 

Let us assume that the chain is very long, so as to neglect boundary effects and consider 
a monomer with variable to in the bulk. The Hamiltonian is then rewritten as 

- PH{i) = fhek ■ z + JbU -Y.^^- PH^'\i) (12) 

iedO 

where —(3H^'^^ (t) stands for the system without the monomer to and dO denotes the set 
of sites neighbouring the site 0. The constrained partition function Z (to) ^'^ that site 
in the bulk reads 

Z [k] = e^'^'"-' [ dHaoe'^'^-^^^^oUQiioo) (13) 

with tgQ = {ti\i G do} and where we have defined 

Q(4o) = / d2t\(ou9o) e-'^^*'"^*) , t\(ou9o) = {iS ^ U 90} (14) 

as the cavity (i.e. in the absence of site 0) partititon function. The latter obviously 
factorises for open chains, or for very long closed ones, having Qitgo) = YliedoQi^i)- 
This allows us to write the following expressions for the partition function Z (t) in terms 
of its cavity counterpart Q (t) in the bulk, viz: 



Q (i) = e^*^*-^ / d^u e^^*-"Q (n) (15) 

Z (i) = e^'^'-' (^j d'u e-^^*"g (n)^ . (16) 

The expression dH]) for the elongation L{f) simply becomes 

L{f) = CB [ dHPit) [i-z] , (17) 

in which -P (t) = Z{t) /Z and Z = J^a d^t ^ (^) ^ normalisation constant. A 
numerical solution of Q (t) can be found by simple iteration of eq. ( IT5l) . Once Q (t) is 
known we can calculate Z (t) and L{f) by using equations ( !T6|) and (fT7|l . respectively. 
In order to fit the experimental measurements we need to calculate the correlation 
functions appearing in ffTTj) efficiently. In the cavity method this boils down to 
calculating the perturbations of P (t) with respect to Jb and 6^- Indeed, by comparing 
(ITT]) and (HZ]) we note that 

N 



^EE((^-^)(^-M>.= / (18) 

i=i k=i ^ ^ ^ 

f ((<.:■«) (4 =/ d^i^(i.i) (19) 
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where the derivatives of P{t) can be expressed as follows: 



dP (i) 1 dz (i) ^ idz 
^ ^ - ^ ^ - p{t) 



X 



(20) 



dx Z dx ^ ' Z dx 

As Z (£) and Q {u) are related via equation (fT6|) . perturbations of Z{t) correspond to 
perturbations of Q [u), viz. 



dZ (i) 
dZ {i) 



2e 



Qii) 



dQ {i) 



B 



f{i-z)Z{i)+2e-f'^'-'Q (t) 



dQ (i) 



with similar equations for the normalisation constant Z: 



dZ 
dZ 



52 



d2£2e-/bfli-^Q (t) 



dQ (i) 



dJf 



dh 



B 



52 



-/(t.f)Z(t)+2e-/^--Q(t)^ 



B 



Perturbations of Q {u) obey the following equations 



dQ (t) 
dQ ji) 

dhn 



fbst-z / ^JbIu 



52 



[i ■ u) Q {u) + 



dQ {u) 
dJ 



B 



f {i ■ z) Q (t) + ef'^'-' I d\ e-^^*-^^ 



52 



db 



B 



(21) 
(22) 

(23) 
(24) 

(25) 
(26) 



To calculate the correlations one proceeds as follows: the set equations f lT5|25|26p for 
{Q{t), djgQ (t) , dbgQ (t) } is solved numerically. Their solution is then used to calculate 
the perturbations of Z (t) and its normalising constant Z, given by eqs fl2m24p . These 
are subsequently used in (l20l) to obtain the perturbations of P(t) which are finally 
plugged in eqs. fll8lll9p to obtain the correlations. 

Note that is possible to show that one can recover the Ihs of eqs. (I18II19I) from the rhs. To 
do this one only needs to write down an explicit solution of d^Q (t) by formally solving 
its corresponding equation via infinite iteration. Indeed, to illustrate let us consider the 
perturbation with respect to Bb- The solution of equation (126|) can be written as follows 

dQ {U) 



dh 



{ii ■ z) Q (U) 



k=0 

This, in turn, yields: 



5-2 



.e=o 



(27) 

^'^ {ii+k+l ■ ^) Q {ii+k+l) 



dZ {U) 



dh, 



f{U-z)Z{i)+2fJ2 



k=0 



52 



i+l+i 



.e=o 



i+k+l) {ti+k+1 ■ z) 



where Z(ti, tj+i, . . . tj+fc+i) is the constrained partition function of a segment of the chain 
of length k + 1: 

Z{k,k+u ■ ■■ii+k+l) = Q (ii) eSti/fest,+ri+JBELo*".+r^.+.+ig (33) 
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Note that using equation (fTB|) we have that: 

/ d tiZ (ti, . . . ti+k+i) = Z (tj+i, . . . ti+fc+i ) • (29) 

This imphes that the normahsation of tj+i, . . . tj+fc+i) is precisely Z and we can 
write P{ti,ti+i, . . .t 

i+fc+l) Z(ti, ■ ■ - t i+k+i)/Z. These results yield the following 
nice and compact expression for dhgP(t): 

dP (i ) 

= fP (U) [U ■ z - (t; ■ z)\ (30) 

CO „ / k \ 

+ 2/ / I JJ^ d^tj+i+£ I P(ti, tj+i, . . . tj+fc+i) {tj+fc+i ■ 2 — (tj+fc+i ■ ^)} 



fc=0 



which can be used to nicely recover the correlation function appearing on the Ihs of eq. 
[T9l These preceding expressions are however impractical and it is preferable to solve 
the aforementioned equations numerically. 



4.2. Double stranded DNA 

This model is solved as before using the cavity method. Although the expressions in 
this bit more involved, we do not linger in minor details and simply report 

the most revelant equations for the curious reader. As in this case we have two states 
per monomer we obtain equations for the constrained partition function Z„ (t) and its 
cavity counterparts {t)\ 




with cr G {B, S}. The elongation L takes the simple form: 

L (/) = / d^t [i ■ z) [cBPB{i) + hsPsii)] (33) 

where Pa (tj = Z^ (t) jZ with Z = X^o- /sa d^t Z^ (t) a normalisation constant. The 
expressions for the propagation of pertubations in this case are fairly tedious and we 
prefer to report their exact expression somewhere else [23] . 



5. The fitting witii experimental measurements and Monte Carlo 
Simulations 

In this section we summarise the results of our fit to experimental measurements 
provided by the "Small Biosystems Laboratory" in Barcelona, both for single stranded 
and double stranded DNA. To check the validity of our analytical results we have also 
performed Monte Carlo simulations of both Hamiltonians ( 1T1I41) using Heat Bath (24] for 
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Figure 1. Fit of ss-DNA stretching measurements taken at the "Smah Biosystems 
Lab" in Barcelona, 100 experimental points.. Experimental data is reported in 
full squares, theoretical solution in straight line while open circles are Monte Carlo 
simulations. Values of the fit: 6^' = 0.63 nm, Jb = 7.68 pN nm, = 1-3 x lO"'^. We 
take room temperature ksT = 4.04738 pN nm. 



the vector variables (combined with Metropohs for the discrete variables in the second 
model) for a system with = 1000 monomers and open boundary condition. Note 
that we do Monte Carlo simulations using the values of already fitted parameters so 
as to check the correctness of our work. Using Monte Carlo simulation for fitting the 
experimental data is a very time-consuming and impractical task. 

5.1. Single stranded DNA 

We first proceed to fit model ([!]) to a set of experimental measurements of ss-DNA 
stretching. The results are reported in Figure [H where we plot experimental results, the 
fit and Monte Carlo simulations. The experiments were performed in a 10 mM NaCl 
solution on a chain of 3000 bases whose distance is estimated to approximately 0.7 nm 
by X ray diffraction. Because of the reduced range of forces (from about to 40 pN) the 
single monomers' extension is not detectable and therefore we keep the inverse Young 
modulus fixed to zero during the fitting. Leaving free would indeed return a small 
value, thus this choice allows us to reduce the number of fitting parameters. The results 
of the fit are a monomer's length b^^^ ~ 0.63 nm and a coupling constant ~ 7.68 pN 
nm, while for cost function we obtain = 1-3 x 10^'^. This result suggests the model's 
monomer corresponds to the actual DNA base distance, at least for the single stranded 
chain. This is very encouraging, since no discrete model so far could capture the real 
scaling of the molecule. As stated before, in fact, the absence of a bending resistance 
results in a monomer's length longer than the actual one (the Kuhn segment), while the 
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Figure 2. Left figure: Fit of ds-DNA stretching measurements taken at the "Small 
Biosystems Lab" in Barcelona, 100 experimental points. Experimental data is reported 
in full squares, theoretical solution in straight line while open circles are Monte Carlo 
simulations. The force range is narrow to avoid any effect linked to the overstretching 
from above 50 pN and to thermal fluctuation from below 4 pN. Values of the fit: 
b^g^ = 3.22 nm, Y = 4570.38 pN, Jb = 56.26 pN nm, 9 x lO'^. Here we have 

assumed room temperature ksT = 4.04738 pN nm. Right figure: Fit of overstretched 
DNA. Experimental points are plotted in filled squares, theoretical solution in straight 
line and Monte Carlo simulations in open circles. The data shows the typical behavior 
of dsDNA stretched by high forces: the double plateau and the jump around 65 pN 
are easily recognizable. The fit is operated over 100 points, starting from / = 1.5 pN. 
The parameters values returned by the fit are reported in table [H = 4 x 10^'*. 



DPC model yielded a ss-DNA monomer's length much smaller than the base distance 
|15] . This means that the qualitative picture of the model can give meaningful and 
expendible quantitative results if a rigorous analysis of the problem is carried out. This 
last point is also confirmed by the more than good agreement between our analytic 
solution and the Monte Carlo simulation. 

5.2. Double stranded DNA 

Figure E] summarises the results of our fit to experimental measurements for ds-DNA. 
For sake of efficiency, the fitting is done in two stages. Firstly, we fit the first model to to 
experimental measurements with forces below 50 pN, that is, below the overstretching 
transition (left figure in|2]). The agreement is beyond any doubt, finding the monomers' 
length of around = 3.22 nm, the Young modulus to be F = 4570.38 pN and 
Jb = 56.26 pN nm. We then proceed to fit the second model to the set of experimental 
measurements up to range of forces to 80 pN, while keeping the values of fe^"*, Jb and 
Y fixed to the previous fitted values. The results of the fit, together with results from 
Monte Carlo simulations are shown in Figure [2] (right figure) and the parameters are 
reported in Table [1] The result is excellent, proving that this approach can realistically 
approximate the physical system: as already mentioned, we find a native monomers' 
length which is of the order of 10 base pairs (the latter being ~ 3.3 A), which interestingly 
corresponds to the DNA helix period. This is much smaller than what is achieved with 
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^(0) 


hs 


Jb 


Jbs 


Js 


IB 


^BS 


3.22 


5.63 


56.26 


26.35 


0.81 


106.61 


-2.75 



Table 1. Parameters of the fit. For the sake of compactness we write J^a = Ja- The 
parameters are given in nm while Ja„i and are given in pN nm. 

the FJC and somehow better than what is obtained in [15], since the very small value 
(less than a base pair) calculated therein has a difficult physical interpretation. The 
denaturated unit length hs is 1.75 times bigger than the native one, giving back the 
known ratio between S-DNA and B-DNA total length: this value conffims that the helix 
unwinds, even if it doesn't give any insight on the melting into bubbles. Finally we get 
from the values of Jb, Jbs and Js that the native strand is much stiffer (and so more 
correlated as well) than the denaturated portion of the molecule. 

6. Conclusions 

In this paper we have revisited the very simple phenomenological model of Storm and 
Nelson which faithfully reproduces the mechanical properties of dsDNA and where the 
monomer's length represents one helix period of the DNA polymer. The model is treated 
with the cavity method which facilities the fitting by easily allowing to understand 
changes of the fitting parameters as perturbations of the cavity equations. 
We believe that the good analytic and numerical results obtained in this paper 
can be the starting point of many further applications in biophysics and molecular 
biology, particularly in heterogeneous polymers, block copolymers and protein helix-coil 
transitions in general. We also think that more realistic yet solvable models can be drawn 
from this study. For instance, since we uncover that a DPCs monomer corresponds to 
a whole helix period, we could incorporate its internal degrees of freedom in such a way 
to study the single base pair's behaviour. 
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Appendix A. Expansion for small forces 

Looking at the expression ([3]) for //(/), we note tliat tlie dependence on tlie force appears 
explicitly through the definition of and implicitly through the thermal average of 
(ti ■ Sy Expanding the latter in powers of / we write: 

{i. ■ z) = {k ■ z)^^, + f (I {U ■ ^>)^_^ + If {U ■ z))^_^ + O . (A.l) 

Noting that (ti ■ ^)j_o is identically zero (as this is Heisenberg model in one dimension 
at zero external field) and keeping the lowest terms in / we can write: 

L = hffx + 0{f) (A.2) 

with 

1 ^ 

i,k=l 

being the magnetic susceptibility at zero external field. Assuming periodic boundary 
conditions x takes the following form (see, for instance, pT]): 

^ _ 1 Jb + Jb cotanh(JB) - 1 
3 Jb — Jb cotanh( Jb) + 1 

Note that eqs. flA.2p and flA.4p give an explicit expression of L(/, Jb), amenable for 
fitting fe^'' and Jb to experimental measurements but only valid in the linear regime of 
forces. 

A similar analysis can be also done in the region for large forces by expanding in powers 
of Jb instead. This is a pedagogical exercise left to the reader. 



